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HIGHLIGHTS 


•  We  model  rechargeable  batteries  for  use  in  BMS  design. 

•  The  rate  capacity  effect  and  the  state  of  charge  are  explicitly  considered. 

•  The  model  is  based  on  an  electrochemical  approach. 

•  A  reduced-order  model  with  electrochemical  sense  is  presented. 

•  This  approach  is  able  to  interpret  the  most  commonly  used  models  in  the  literature. 
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There  are  numerous  models  of  rechargeable  batteries  in  the  current  literature.  Some  of  them  are 
complicated  electrochemical  approaches;  others,  are  given  by  simple  analogies.  However,  simple  models 
with  electrochemical  states  capable  of  describing  important  quantities,  like  the  rate  capacity  effect  and 
the  recovery  effect,  are  hard  to  find.  In  this  paper,  based  on  an  electrochemical  approach,  we  present  a 
generalized  model  suitable  for  use  in  BMS  applications,  which  takes  into  account  explicitly  the  rate 
capacity  effect  and  the  state  of  charge.  Moreover,  the  model  is  thought  up  for  general  energy  storage 
processes  based  on  mass  transport  and  charge  transfer.  The  proposed  general  model  approach  is  able  to 
interpret  the  most  commonly  used  models  in  the  literature. 

©  2013  Elsevier  B.V.  All  rights  reserved. 


1.  Introduction 

The  battery  model  is  necessary  both  for  estimation  and  pre¬ 
diction,  which  are  very  important  topics  for  battery  management 
systems  (BMS).  A  good  model  is  needed  in  order  to  know  the  on¬ 
line  state  of  charge  (SOC),  the  state  of  health  (SOH),  and  the  avail¬ 
able  power,  also  to  decide  control  actions  for  equalization  and  for 
battery  pack  design.  An  important  feature  of  rechargeable  batteries 
is  that  the  voltage  drops  faster  for  high  discharge  currents.  Thereby, 
the  effective  capacity  drops  for  high  discharge  rates.  This  effect  is 
called  rate  capacity  effect  (RCE)  and  it  is  associated  with  the  re¬ 
covery  effect  (RE),  which  is  the  time  required  for  battery  recuper¬ 
ation.  Both  quantities  are  important  for  an  optimal  use  of  batteries. 
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In  the  current  literature  there  are  numerous  models  of  rechargeable 
batteries.  Some  of  them  are  complicated  electrochemical  ap¬ 
proaches;  others,  are  given  by  simple  analogies.  However,  it  is 
difficult  to  find  practical  models  that  interpret  the  RCE  and  the  RE 
using  electrochemical  states.  In  this  paper,  based  on  an  electro¬ 
chemical  approach,  we  present  a  generalized  model,  suitable  for 
use  in  BMS  applications,  which  takes  into  account  explicitly  the  RCE 
and  the  RE.  The  novel  contribution  is  that  our  model  has  a  simple 
structure  composed  by  the  cascaded  of  a  pure  integrator,  a  high 
pass-filter  (HPF),  and  a  nonlinear  function;  all  these  components 
related  with  explicit  important  electrochemical  states,  like  SOC, 
and  surface  charge  concentration,  which  is  a  key  variable  for 
describing  the  concepts  of  RCE  and  RE.  Moreover,  the  model  is 
proposed  for  general  energy  storage  processes  based  on  mass 
transport  and  charge  transfer.  This  approach  is  able  to  interpret  the 
most  commonly  used  models  in  the  literature. 

The  development  of  an  accurate  dynamic  battery  model  is 
difficult  due  to  the  nonlinear  infinite-dimensional  and  distributed 
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mass  transport  processes  governing  electrochemical  system  dy¬ 
namics.  Physiochemical  models  provide  an  accurate  insight  on 
battery  processes,  but  models  based  on  fundamental  physi¬ 
ochemical  principles  are  rarely  employed  in  real  time,  due  to  their 
high  order  and  computational  requirements.  The  need  to  identify  a 
great  number  of  parameters  is  also  a  drawback.  Several  kinds  of 
simplifications  have  been  presented  in  some  publications;  for 
example,  diffusion  limitations  in  the  electrolyte  and  potential  dis¬ 
tribution  in  the  porous  structure  of  the  electrodes  are  disregarded 
in  Ref.  [1],  Other  models  are  focused  on  the  simplification  of  the 
diffusion  dynamics  in  the  active  material.  In  this  direction,  in  Ref. 
[2]  a  second  degree  polynomial  approximation  is  used  to  model  the 
concentration  profile  in  the  active  material.  Barbarisi  [3]  presents  a 
Ni-MH  electrochemical  model  driven  by  the  measured  current  and 
potential.  For  this  model  several  kinetic  and  transport  parameters 
of  the  battery,  which  are  difficult  to  obtain,  need  to  be  known.  In 
order  to  overcome  this  difficulty,  the  static  nonlinearity  can  be 
approximated  by  using  a  Taylor  series  expansion  and  by  identifying 
the  coefficient  applying  linear  regression,  as  in  Ref.  [4],  In  recent 
publications  [5—7]  reduced  order  models  based  on  frequency 
domain  transfer  function  analysis,  are  presented. 

A  widely  used  approach  to  model  batteries  consists  in 
employing  analogies  with  electrical  circuits  where  the  voltage 
source  is  the  open  circuit  voltage  (OCV).  The  analogies  with  elec¬ 
trical  circuits  are  presented  in  Refs.  [8,9]  as  well  as  the  references 
they  contain.  These  models  have  been  extensively  used,  for 
example  in  Ref.  [10],  to  estimate  the  SOC  in  a  Lithium  Polymer 
battery  with  nonlinear  sliding  observers;  in  Ref.  [11],  for  hybrid 
electric  vehicles  with  lead  acid  batteries;  in  Ref.  [12],  with 
nonlinear  observers  to  estimate  the  parameters  of  the  model  and 
the  SOC  jointly;  also  in  Ref.  [13],  to  estimate  the  SOC  in  lead  acid 
batteries.  Such  models  are  simple,  but  difficult  to  relate  with  the 
battery  chemical  states.  One  well-known  model  used  to  represent 
the  dynamics  of  the  stored  charge  for  design  engineering  applica¬ 
tions  was  presented  in  Ref.  [14],  It  is  a  battery  model  based  on  the 
approach  of  chemical  kinetics,  known  as  kinetic  battery  model. 
Specifically,  it  takes  into  account  the  rate  capacity  effect  by 
assuming  that  the  charge  can  be  stored  in  two  ways,  either  as 
immediately  available  or  as  chemically  bound.  Although  this  model 
interprets  the  RCE,  is  not  able  to  represent  the  output  voltage  as  a 
natural  function  of  electrochemical  states  requiring  special  ad¬ 
justments  based  on  tables. 

Basically,  simple  analogies  are  practical  but  empirical  and 
sometimes  they  do  not  work  properly.  On  the  other  extreme,  there 
are  electrochemical  models  with  several  parameters  and  states, 
which  are  very  difficult  to  obtain.  In  this  paper,  following  electro¬ 
chemical  arguments,  we  show  that  any  rechargeable  battery  can  be 
modeled  as  a  cascade  of  a  linear  dynamic  system  with  a  static 
nonlinearity  in  a  setup  called  Wiener  model.  The  linear  dynamic 
system  is  related  to  the  transport  of  the  reacting  substances  in  the 
active  material.  These  processes  are  governed  by  Fick’s  law  and  can 
be  represented  as  a  series  of  two  subsystems;  one  is  a  pure  inte¬ 
grator  representing  the  real  stored  charge  or  SOC,  and  the  other  is  a 
high-pass  filter  representing  the  dynamic  of  the  concentration  at 
the  electrode  interface  and  related  to  the  rate  capacity  of  the  bat¬ 
tery.  The  static  nonlinearity  is  due  to  the  electrochemical  reactions 
at  the  electrode  interfaces,  and  it  is  governed  by  a  Butler-Volmer 
type  equation  but  can  also  be  combined  with  other  reaction  types 
described  by  a  non-linear  function  which  characterizes  the  elec¬ 
tromotive  force  of  the  battery. 

The  paper  is  organized  as  follows:  in  Section  2  we  present  the 
model  formulation  basically  through  two  main  stages.  First,  the 
charge  transfer  process  at  the  interface;  and  second,  the  diffusional 
process.  In  this  section,  the  diffusional  process  is  shown  to  be 
composed  by  two  parts:  one,  a  pure  integrator,  and  the  other  a  HPF. 


The  general  model  is  presented  in  this  section.  In  Section  3,  the  rate 
capacity  effect  is  related  to  the  diffusional  process.  In  Section  4,  we 
present  a  reduced  order  approximation  and  it  is  shown  that  in  the 
case  of  second  order  approximation,  the  dynamics  of  diffusion 
follow  the  kinetic  battery  model.  Connections  with  the  Thevenin 
electrical  circuit  analogy  are  also  presented.  In  Section  5,  a  proce¬ 
dure  to  identify  the  parameters  together  with  experimental  results 
are  shown  for  different  battery  types,  like  Li-Ion,  Ni-MH,  and  Lead 
Acid.  Finally,  we  discuss  the  conclusions  in  Section  6. 

2.  Model  formulation 

A  battery  is  basically  a  system  conformed  by  two  electrodes 
immersed  in  an  electrolytic  media.  During  discharge,  the  negative 
electrode  is  oxidized,  while  the  positive  electrode  is  reduced.  The 
process  is  reversed  during  charge,  being  this  global  process 
responsible  for  energy  storage  or  release.  Battery  electrodes  are 
generally  porous  structures  composed  of  particles  of  active  mate¬ 
rial,  conductive  additives  and  adequate  binders,  supported  on  a 
conductive  substrate  acting  as  current  collector,  the  whole  assem¬ 
bly  being  immersed  in  the  electrolytic  media.  Models  for  porous 
electrodes  mainly  use  the  porous  electrode  theory,  as  reviewed  by 
Newman  and  Tiedemann  [17],  In  according  to  this  theory,  the  liquid 
and  the  solid  phases  are  approximated  as  superimposed  continua. 
The  intimate  contact  between  both  phases  generates  the  interfacial 
surface,  where  the  different  electrochemical  processes  take  place 
(electrochemical  interface).  In  the  energy  storage  or  release  pro¬ 
cesses,  two  main  stages  must  be  distinguished:  one,  corresponding 
to  the  charge  transfer  process  at  the  electrodes  electrochemical 
interface,  including  electrochemical  reactions  and  double  layer 
charging;  and  another  one,  corresponding  to  mass  transfer,  either 
in  the  electrolyte  or  in  the  electrode  active  material.  In  the  next 
subsections,  we  describe  both  stages,  which  are  the  main  part  of 
our  battery  model. 

2.1.  Faradaic  charge  transfer  processes 

The  electrochemical  reactions  taking  place  at  the  electro¬ 
chemical  interface  give  rise  to  the  faradaic  current,  I(t).  The  faradaic 
current  may  be  described  in  terms  of  the  Butler-Volmer  equation, 
as  follows  [18,19]: 

i = kV£(nc?) + krear£(  nc0>  (i) 

where  k  and  k  are  constants  corresponding  to  reduction  and 
oxidation  reactions;  E  is  the  electrode  potential;  ar  —  arFIRT  and 
a°  =  a°FIRT  are  constant  with  ar  =  a°  -  1  and  a0  is  a  symmetry 
factor  in  the  interval  (0,1 );  F  is  the  Faraday  constant;  R  the  universal 
gas  constant;  T  is  the  temperature;  and  c-  and  c°  are  the  concen¬ 
tration  of  reactants  and  products  at  the  electrochemical  interface. 
In  the  sequel  we  shall  consider  that  the  very  high  concentrations  of 
species  are  constants,  so  for  most  batteries  it  is  possible  to  reduce 
the  complexity  to  a  simple  redox  reaction,  as  follows  [18,20]: 

I  =  k°c°e“°£  -  ii  crea'E.  (2) 

By  multiplying  and  dividing  the  first  term  on  the  right  by  c°,  and 
by  cr  the  second,  which  are  the  maximum  concentration  of  the 
oxidized  and  reduced  states  of  the  active  materials,  the  following 
simple  equation  holds: 

I  =  k°X°e“°£  -  krXrea'E,  (3) 

where  kr  =  k  ?  and  k°  =  k  c°  are  constants,  and  Xr  —  cr  jd  and 
X°  =  c°/c°  are  the  normalized  concentrations  at  the  interface,  with 
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respect  to  the  maximum  concentration.  Taking  into  account  that 
concentration  lies  in  the  interval  [0-1],  it  follows  that  each  elec¬ 
trode  fulfills  Xr  =  1  -  X°.  In  what  follows,  for  notation  simplicity,  we 
neglect  the  superfix  o  by  calling  X  =  X°,  c  =  c°,  and  c  =  c°.  Then, 
equation  (3)  can  be  written  as 

I  =  fi{E,X)  =  k°Xea°E  -  kr(  1  -X)ea'E,  (4) 

where  fi(EJ()  is  the  Butler— Volmer  equation  expressing  the  func¬ 
tional  relationship  of  current  with  respect  to  potential  and  con¬ 
centration.  From  equation  (4)  it  is  possible  to  obtain  the 
concentration  at  the  electrode  surface  X  explicitly,  as 

uearE  ,  i 

■X  =  fx(E,  I)  =  for  gars  +  fcogg°E  •  (5) 

Note  that  concentration  X  increases  when  the  current  that  enter 
to  the  electrode  is  considered  positive.  Although  it  is  not  possible  to 
obtain  an  explicit  expression  of  the  form  E  =  from  equation 
(5),  in  the  Appendix  A  it  is  demonstrated  that  such  function  exists. 
It  is  important  to  note  that  when  the  faradaic  current  is  zero,  (/  =  0), 
by  using  logarithms  in  equation  (5),  the  potential  can  be  obtained 
explicitly  as 

<6> 


with  the  electrolyte. 


N  small  volumes,  each  with  volume  (V)  equal  to  each  other,  with 

constant  fractional  concentrations  denoted  as  x,-  for  i  =  1 . N.  Then, 

by  denoting  the  mass  flow  as  Jy,  in  mass/time  units,  from  com¬ 
partments  i  to  j,  the  following  relationships,  Vj  =£  i,  hold: 
lij  =  Kjj(Xi  -  Xj),  (8) 


which  is  called  Electro-Motive-Force  (EMF).  The  battery  EMF  is 
denoted  as  E  — /e(0,X),  which  represents  the  functional  relationship 
of  the  potential  with  respect  to  the  concentration  X  at  zero  current, 
1=0.  Note  that  the  potential  £  increases  when  the  concentration  X 
increases. 

It  is  worth  noting  that  the  potential  at  the  electrode  surface,  E(t), 
can  be  decomposed  into  a  sum  of  two  terms.  One,  is  the  equilibrium 
potential,  defined  by  zero  current,  I  =  0,  and  concentration 
Xeq  =  0.5,  denoted  as  Eo  and  is  obtained  using  equation  (6).  The 
other,  is  the  potential  that  depends  on  both  the  current  and  the 
concentration. 

2.2.  The  diffusion  process 

The  reacting  substance  approximates  the  electrode  surface  by 
diffusion.  The  second  part  of  the  model  deals  with  the  dynamics  of 
the  mass  transport  process,  which  determines  the  concentration 
X(t)  described  in  equation  (5).  These  are  essentially  diffusion  pro¬ 
cesses,  either  in  the  electrolyte  or  in  the  electrode  active  materials. 
The  number  of  substance  moles  that  reacts  at  the  electrochemical 
interface  is  proportional  to  the  electrons  produced  or  consumed, 
according  to  the  Faraday  constant.  The  electrons  flow  through  the 
load  while  the  reacting  substance  is  accumulated.  The  total  amount 
of  such  moles  is  proportional  to  the  accumulated  amount  of  charge 
giving  the  SOC  of  the  battery,  as  follows: 


where  Qo  and  Qmax  are  the  initial  and  the  maximum  battery  charge. 
The  quantity  Qmax  is  called  battery  capacity  and  is  defined  as  the 
product  of  a  very  low  constant  current  multiplied  by  the  time 
required  to  reach  the  maximum  charge,  starting  with  the  battery 
completely  discharged  (Qo  =  0). 

In  order  to  model  the  diffusional  process,  consider  that  the 
reacting  substance  is  accumulated  or  released  by  diffusion  from  an 
area  of  higher  concentration  to  another  of  lower  one.  Assume  that 
the  physical  place  where  the  reacting  mass  is  stored  is  divided  into 


Kji  =  Ky  (9) 

Jy  =  - Jji ;  (10) 

where  Ky  is  in  mass/time  units;  suffix  1  =  0  corresponds  to  the 
differential  compartments  adjacent  to  the  contact  surface  with  the 
electrolyte,  the  interface.  There  is  also  an  associated  mass  balance, 
which  is  given  by 

with  y  =  l /Vc  in  1/mass  units,  and 

E/i-p.  (12) 

where  M  is  the  number  of  compartments  adjacent  to  the  surface  in 
contact  with  the  electrolyte.  The  state  of  charge  is  the  total  accu¬ 
mulated  charge  and  it  is  given  by 

soc(t)  =  Ensr-  («) 

i=i, jv  ,v 

A  simple  example  of  three  compartments  will  be  shown  for 
clarity.  Assume  that  the  material,  where  the  mass  is  stored,  has 
N  =  3  with  M  =  2  compartments,  as  it  is  shown  in  Fig.  1.  Then,  the 
mass  balance  equations  are  given  by 

=  y(Joi  ~ J 12  —  J13); 

* 2  =  y(Ju+J32), 

*3  =  y(Jo3  +Jl3  -J32), 

where 


J12  =  ^12(Xl  -x2), 

Jl3  =  Kl3(Xl— x3),  (14) 

J32  =  1^32  (X3  —  x2). 
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Joi  and  J03  are  inputs.  Taking  into  account  that  Ky  —  Ky,  and  defining 
ai  and  d3  as  the  fractional  part  of  the  sum  of  incoming  mass  flow 
/  =  Jm  +  J03.  such  that  a,  +  a3  =»  %  the  following  holds: 

Joi  =  ai/;  Jo3  =  a3JC-  (15) 

The  system  equations  can  be  written  in  a  matrix  form  by 

x(t)  =  Ax(t)+yBf,  (16) 

where 


Fig.  2.  Battery  model. 


X(t)  =  [xi,x2,x3]r;  B  =  [a3 , 0, a3]r 

/-K12-/C13  Ku  K13  \ 

A  =  r  Ku  -K32-Ku  K32 

\  K13  K32  -K32  -  Ku ) 

We  are  interested  in  obtaining  the  dynamical  relationship  be¬ 
tween  the  charge  current  and  the  concentration  at  the  interface  X. 
To  this  end,  let  us  average  the  individual  concentrations  corre¬ 
sponding  to  the  M  differential  volumes,  at  the  interface,  in  contact 
with  the  electrode  surface.  This  is  obtained  by  multiplying  BT  by  the 
vector  x,  as  follows: 

X  =  Brx(t).  (17) 

The  transfer  function  between  concentration  and  charge  current 
is  obtained  from  this  definitions  by  applying  the  Laplace  transform 
to  equations  (16)  and  (17)  resulting  in  the  following  input/output 
relationship: 

X(S)  =  H(S)I(S),  (18) 

where  X(S)  and  1(S)  are  the  Laplace  transform  of  X  and  /;  S  is  a 
complex  variable,  and 

H(S)  =  BT(Sl-Ay1By/F,  (19) 

where  I  means  identity  matrix.  Before  continuing,  it  is  important  to 
say  that  the  same  procedure  can  be  extended  to  the  case  of  N 
constant  volume  partitions,  as  follows: 


x(f)  = 

-  ;[x1,x2,---,xjv]r 

;  /  = 

=  phi'. 

/-.£  Ky 

Ku 

'”V13  ... 

Kw  \ 

A  = 

]~k2, 

.E 

I<2]  k23  ... 

K2n 

Km 

Kn 2 

...  kN3  - 

E  KNj 

\ 

i=t,JV 

where 

Joi  represents  the  mi 

iss  flow  that  enter 

to  each  ol 

compartments  in  contact  with  the  surface,  and  B  is  a  column  vector 
having  a  constant  entry  value  a,  if  the  infinitesimal  volume  V  is  in 
contact  with  the  surface  and  zero  otherwise.  The  M  constants  a,  are 
constrained  to  J2  a,-  =  1.  The  total  mass  flow  generated  by  the 
i = 1  ,M 

reacting  substances  at  the  surface  of  the  electrode  is 
/  =  J2  Joi  =  J2  ajc.  Note  that  ideally,  the  expanded  matrices 

i  —  1  ,M  i  — 1,M 

A,  B,  and  C  are  infinite-dimensional,  but  still  can  be  used  in  equa¬ 
tions  (17)— (19),  obtaining  an  infinite-dimensional-linear  dynamic 
system. 

By  multiplying  and  dividing  equation  (18)  by  SQmax,  the  func¬ 
tional  relationship  between  the  current  and  the  concentration  at 
the  interface  can  be  rewritten  as  a  function  of  the  state  of  charge  as 


X(S)  =  G(S)SOC(S), 

(21) 

where 

G(S)  =  H(S)SQmax; 

(22) 

(23) 

2.3.  The  complete  model 

It  should  be  noted  that  the  battery  is  formed  by  two  electrode¬ 
electrolyte  interfaces.  However,  since  generally  only  one  is  limiting, 
it  is  possible  to  represent  the  global  dynamics  with  only  one 
equivalent  electrode.  Apart  from  diffusion  and  static  nonlinearity, 
the  complete  battery  model  should  include  the  electrolyte  ohmic 
resistance  (re)  and  the  effect  of  the  double  layer  capacity  (Cdi).  The 
complete  model  is  obtained  by  relating  the  charge  transfer 
E  =  /e(/,X),  the  diffusion  processes  (21 ),  re,  and  the  current  (I a), 
which  includes  the  current  due  to  the  capacity  of  double  layer  and 
others  derived  from  secondary  reactions  such  that  /bat  =  1  +  1 a, 
where  /bat  is  the  current  measured  at  the  battery  terminals.  Thus, 
the  rechargeable  battery  can  be  modeled  as  a  cascade  of  a  linear 
dynamic  system  followed  by  a  static  nonlinearity,  as  shown  in 
Fig.  2,  which  is  known  as  the  Wiener  model  [4], 

3.  The  rate  capacity  effect 

Since  the  redox  reaction  occurs  at  the  electrode  surface  and  that 
the  charge  transport  process  is  diffusive,  for  a  given  charge/ 
discharge  current,  X  increases/decreases  faster  than  the  SOC.  By 
taking  into  account  that  E  depends  directly  on  X,  the  useful  capacity 
of  the  battery  heavily  depends  on  the  charge/discharge  current 
profile.  Great  current  demands  reduce  the  effective  capacity  of  the 
battery.  This  effect  is  reflected  on  the  high-rate  discharge  behavior 
of  the  battery  and  it  is  called  rate  capacity  effect  (RCE),  [15],  and  it  is 
important  for  battery  use  optimization,  [16],  It  can  be  seen  from 
equation  (21)  that  X  and  SOC  are  related  by  the  dynamics  imposed 
by  the  filter  G(S),  which  determines  the  rate  capacity  H(S)of  the 
battery  and  also  the  recovery  time.  Thus,  we  are  interested  in 
studying  the  filtering  properties  of  G(S). 

First,  note  that  G(S)  is  an  infinite  dimensional  rational  function 
in  the  variable  S.  Secondly,  its  steady  state  gain  is  equal  to  one.  To 
prove  the  last,  note  that  in  steady  state  the  mass  is  equally 
distributed  in  all  the  small  volumes,  and  by  considering  that  both 
the  SOC  and  the  fractional  X  are  in  the  interval  [0,1],  it  follows  that 
in  steady  state  both  quantities  are  equal.  Thus,  from  equation  (21), 
it  means  that  the  steady  state  gain  of  G(S)  is  equal  to  one,  G(0)  =  1. 
Now  we  will  show  that  GQ'w)  is  an  infinite  dimensional  rational 
HPF.  To  this  end,  let  us  define  a  positive  real  function  and  two 
important  properties. 
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A  rational  function  H(S)  is  called  positive  real  function  —  PRF  —  if 
and  only  if  H(S)  is  real  when  S  is  real,  and  Re(H(S))  >  0  when 
Re(S )  >  0,  see  Ref.  [21]  for  details.  If  H(S)  is  a  PRF,  it  is  called 
impedance  and  the  following  properties  hold:  i)  the  absolute  value 
of  the  phase  angle  of  H(S)  is  never  greater  than  w/2;  ii )  H(S)  cannot 
have  any  poles  or  zeros  in  the  right-half  S-plane,  i.e.  it  is  a  minimum 
phase  system.  In  the  Appendix  B  it  is  shown  that  H(S)  in  equation 
(18)  is  a  PRF.  Then,  because  H(S )  is  a  cascade  of  a  pure  integrator 
element  and  the  filter  G(S),  it  is  clear  that  in  order  to  fulfill  the 
property  i,  the  element  G(S)  must  have  a  positive  phase.  From  the 
property  ii,  and  taking  into  account  the  foregoing,  the  module  must 
necessarily  be  greater  than,  or  equal  to,  unity.  Given  that  low- 
frequency  gain  is  unity,  the  filter  frequency  response  of  G(jw ) 
necessarily  amplifies  high  frequencies.  Thus,  G(/w)  is  a  HPF. 

As  an  example,  consider  the  case  where  a  planar  diffusion 
process  is  governed  by  Fields  law,  which  is  a  very  good  approxi¬ 
mation  for  some  batteries.  The  Laplace  transform  of  concentration 
X(S)  is  derived  in  Ref.  [22],  In  the  case  of  finite  diffusion  and 
impermeable  barrier,  it  has  the  following  expression: 


X(S) 


L  coth  </>(S) 
FAD  <p>(SJ 


X(S), 


(24) 


where  D  is  the  diffusion  constant,  L  and  A  are  the  thickness  depth 
and  the  area  of  the  electrode,  and  <p(S)  =  L^JS/D.  Multiplying  and 
dividing  by  <p(S )  the  expression  can  be  rewritten  as 


*  =  (SG(S)’  (25) 

where  Qmax  =  1 IFV  with  volume  V  =AL  and 


G(S)  =  <pCOth  cp,  (26) 

The  frequency  response  G(S  =  jw)  corresponds  to  a  linear  infinite 
dimensional  HPF,  unity  static  gain,  with  frequency  response  as 
shown  in  Fig.  3.  A  similar  analysis  can  be  made  for  spherical 
diffusional  processes.  In  such  cases,  the  concentration  has  the 
following  Laplace  transform: 


X(S) 


L  tanh  q>(S) 

FAD  (tanh  <p(S)  -  <p{S))  (  h 


(27) 


KD2: 


co  L2/D 


where  L  is  now  the  radio  and  A,  the  area  of  the  sphere.  By  consid¬ 
ering  that  the  volume  of  the  sphere  is  V=AL/ 3,  the  equation  can  be 
written  as  equation  (25),  where  G(S)  is  now  given  by 


<p2(S)  tanh  ,p(S) 

3  (tanh<p  (S)  -  <p(S))’ 


(28) 


with  frequency  response  also  shown  in  Fig.  3.  Note  that,  in  both 
cases,  the  function  admits  infinite  derivatives  with  respect  to  S. 
Thus,  they  can  be  written  as  an  infinite  dimensional  rational  HPF  in 
S,  as  in  equations  (19),  (20)  and  (22). 


4.  A  practical  approximated  model 


As  it  was  shown,  the  filter  G(S)  is  an  infinite  dimensional  linear 
system  described  by  ordinary  differential  equations.  However,  in 
order  to  have  a  generalized  working  expression,  we  assume  that  a 
finite  order  can  be  used  as  a  good  approximation.  The  following 
unit  gain  nth-order  rational  function  is  used  to  approximate  the 
infinite  dimensional  G(S): 


C(S)  = 


N(S)  +  n„S" 
D(S)  ’ 


where 


(29) 


N(S)  =  d0  +  niS  H - h 

D(S)  =  d0  +  d]S  H - h  d„_iSn_1  +  S". 

Note  that  the  HPF  requires  equal  order  for  numerator  and  de¬ 
nominator  polynomials  and  equal  constant  term  do,  for  unit  gain. 
By  replacing  the  above  in  equation  (21),  the  concentration  can  be 
written  conveniently,  for  further  parameter  identification  pur¬ 
poses,  as  the  sum  of  two  terms,  one  depending  on  SOC  and  the 
other,  depending  on  I,  as  follows: 


X(S)  = 


N(S) 

D(S) 


SOC  (S)  + 


n„Sn-i 

QmaxD(S) 


X(S). 


(30) 


Any  given  transfer  function  which  is  strictly  proper  (denomi¬ 
nator  polynomial  order  greater  than  numerator  polynomial  order), 
can  easily  be  transferred  into  state-space  observable  canonical  form 
[32],  Then,  a  useful  way  to  write  the  above  equation,  to  be  used  in 
what  follows,  is  in  the  time-domain  as 

x(t)  -  $x(t)  +  r1soc(t)  +  r 2/(t),  (3i) 

where 


(-dm  1 

0 

0  ••• 

o\ 

0 

1 

0 

0 

1\ 

0 

:  1'.  : 

-di  0 

0 

0 

1 

{  d0  J 

\  -do  0 

0 

0  ••• 

o) 

and  the  concentration  X(t)  is  given  by 

m  =  Cx,  (32) 


Fig.  3.  Frequency  response  of  filter  C(j !io)  for  planar  diffusion  symmetry  in  solid  where  C—  [1,0,.  ..0].  Using  the  above,  the  Complete  model  in  State 
electrode.  variables  can  be  written  as  follows: 
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Fig.  4.  Model  approximation  where  HPF  is  the  finite  order  unit  gain  high  pass  filter 
approximation  of  G(S)  representing  the  rate  capacity  effect. 


Re 


Fig.  5.  Equivalent  circuit  for  second  order  approximation.  C  =  (W/p,  Re  =  p(n,  - 1)/ 
doQmax  and  Ce  =  Qmjpfa  - 1). 


(33) 


Ebat(t)  =/E(/;Cx)  +  re/bat(t),  (34) 

where  x  is  a  disturbance  vector  to  the  states  which  includes  the 
effect  of  current  I\  and  other  possible  disturbances.  £|Mt  is  the  po¬ 
tential  at  the  battery  terminals  which  include  the  electrode  po¬ 
tential,  E,  plus  the  ohmic  drop. 

There  are  several  approaches  that  can  be  used  to  choose  the 
appropriate  model  order  (n).  Statistical  methods,  such  as  the  Akaike 
criterion,  [33],  or  the  robust  approach  used  in  Ref.  [6],  The  reduced 
order  approximation  of  HPF  together  with  a  normalized  basis  for 
the  nonlinear  voltage  equation,  like  polynomials,  expansion  of 
Taylor  series,  wavelets,  radial  basis  functions,  etc.,  could  be  used  for 
identification  of  the  Wiener  model  (33)  and  (34)  (see  Refs.  [25,26]). 
However,  a  further  approximation  can  be  made  by  considering  the 
voltage  model  as  formed  by  the  sum  of  two  terms.  One  nonlinear 
that  depends  solely  on  the  concentration  and  another  linear,  that 
depends  only  on  the  current: 

Ebat  =  fE(I,X)  +  re/bat  »/«('  0,X)  +  req/bat,  (35) 

where  the  equivalent  resistance,  req  =  dEbat/9/bat.  includes  not  only 
the  resistance  of  the  electrolyte,  but  also  the  resistance  of  the 
charge  transfer.  In  Appendix  C  a  formal  expression  of  req  is  given. 
Although  the  resistance  due  to  charge  transfer  is  dependent  on 
both  current  and  state  of  charge,  it  is  considered  constant. 

It  is  important  at  this  point  to  establish  the  relationship  between 
EMF  (/e(0,X))  and  the  OCV  (f£( O.SOC)).  The  EMF  is  the  function 
E  =  /g(0,X)  which  relates  the  potential  E  with  concentration  X  and 
the  OCV  is  the  same  as  EMF  but  considered  in  steady  state.  Since  at 
steady  state  fulfills  X  =  SOC,  in  such  conditions,  both  EMF  and  OCV 
are  the  same.  The  final  approximated  model  is  depicted  in  Fig.  4, 
where  the  disturbance  vector  to  the  states  x  was  neglected, 
consequently  I  =  fbat  fulfills. 


4.3.  The  first  order  HPF  approximation  and  its  relationship  with 
other  known  linear  models 


Considering  a  first  order  approximation  of  the  high-pass  filter  in 
equation  (33)  and  neglecting  the  disturbance  vector  x,  we  obtain 
the  following  second  order  model  for  diffusion: 


In  order  to  obtain  a  linear  expression  for  the  voltage  equation,  a 
further  approximation  can  be  made  from  equation  (34)  by  linear¬ 
izing  the  increments  of  the  nonlinear  term,  as  follows: 


Ebat  =  fE(  0,X)  +  req/bat  =  E0  +  V*  +  What,  (37) 


where  p  and  Eo  are  constants  used  for  the  pseudolinear  approxi¬ 
mation  of  the  EMF.  By  using  this  pseudolinear  model,  it  is  possible 
to  derive  an  analog  electric  circuit  by  performing  the  Laplace 
transform  to  equations  (36)  and  (37),  which  gives 


ws>-10 


+  req/bat(S) 


Ebat(S) 


\SQmax 


j.  $"1  -1) 

(dg  +  S)Qmax 


+  Ec 


(38) 

(39) 


This  expression  can  be  represented  by  the  circuit  of  Fig.  5.  Note 
that  the  capacitance  is  related  to  the  maximum  battery  charge  by 
C  =  Qma x/m  and  the  parallel  RgCg  represents  the  first  order  HPF 
approximation,  where  Rg  —  p[n\  -  l)/doQmax  and  Cg  —  Qma x/ 
p{n-[  -  1 ).  Often,  the  series  of  E0  and  the  capacity,  C,  is  replaced  by 
the  OCV.  The  circuit  structure  is  called  Thevenin  Battery  Model  and 
can  be  extensively  found  in  the  literature  (see  for  example  Refs. 
[29,30]  used  for  State  of  Charge  estimation;  also  computational 
implementation  costs  with  respect  to  higher  order  degree 
approximation  were  studied  in  Ref.  [31]).  It  is  noteworthy  that  the 
circuit  is  analog  to  our  model,  from  battery  current/potential  point 
of  view,  only  if  the  OCV  is  approximated  by  the  series  of  a  capaci¬ 
tance  and  a  constant  voltage  source,  which  constitute  a  quasi-linear 
system,  otherwise  the  nonlinearity  of  the  OCV  affects  different  in 
both  models. 

Other  popular  battery  model  used  to  represent  the  dynamics  of 
the  stored  charge  for  design  engineering  applications  is  presented 
in  Ref.  [14],  It  is  a  battery  model  based  on  the  approach  of  chemical 
kinetics,  known  as  kinetic  battery  model  also  known  as  the  hydraulic 
model;  and  it  is  depicted  in  Fig.  6.  Specifically,  it  takes  into  account 
the  rate  capacity  effect  by  assuming  that  the  charge  can  be  stored  in 
two  ways,  either  as  immediately  available  or  as  chemically  bound.  It 
is  generally  described  by  the  analogy  of  two-tanks.  The  bound 
charge  can  be  released  at  a  rate  proportional  to  a  constant  g.  The 
total  charge  (Qmax)  in  the  battery  at  any  time  is  the  sum  of  both  the 
immediately  available  (qi),  and  the  chemically  bound  (<72)- 
Considering  that  q\  and  q2  are  the  tank  levels  of  areas  and  1  -  /? 
respectively,  the  equations  describing  the  charge  dynamics  are 
given  by 


Fig.  6.  Equivalent  hydraulic  model  or  kinetic  battery  model. 
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mt)  =  (q2(t)  -  qi(t))g,  (40) 

(1  -/J)q2(t)  =  (<h(t)  -  q2(t))g  +  u(t),  (41) 

where  /?  e  [0,1]  is  the  ratio  of  available  charge  capacity  to  total 
capacity  and  g  is  a  rate  constant  in  1/time.  By  performing  the 
following  change  of  coordinates: 

SOC  =  0qt  +  (1  -  U)q2,  (42) 

X  =  q2,  (43) 

the  following  equations  succeed: 


(44) 

which  is  equal  to  equation  (36)  with  the  equivalences  do  =  g/ 
/?(1  -  n i  =1/(1  -  /?),  and  u  =  /bat/Qmax-  Note  that  ni  is  lesser  that 
one,  as  it  is  expected  according  to  the  characteristic  of  the  HPF.  Note 
that  this  dynamic  is  equal  to  the  linear  second  order  approximation 
(36)  of  our  general  model.  This  dynamical  part  of  the  model  needs 
to  be  completed  with  the  battery  voltage  model.  A  voltage  model 
given  by  £bat  =  E0  +  kip  +  k2pl(k3  +  p)  +  reqlbat  was  proposed, 
where  klf  k2,  k3  are  constants,  and  p  is  the  normalized  charge 
removed  from  the  battery  that  is  obtained  from  the  discharge  data, 
as  explained  in  Ref.  [27],  Sometimes,  other  voltage  models  are  used, 
such  that  the  Unnewehr  Universal  Model  [28]  expressed  by 
Ebat  =  Eo  —  kSOC  +  req/bat,  where  k  is  a  constant.  However,  in  all 
these  cases  the  voltage  model  is  different  from  equation  (34),  since 


we  used  X  instead  of  SOC,  which  is  electrochemically  correct.  Thus, 
with  this  model  approach,  the  RCE  model  corresponds  with  our 
second  order  model  approximation,  but  the  voltage  model  is 
built  using  tables  which  is  not  practical  for  BMS  applications. 
Others  voltages  models  uses  SOC  instead  of  X  which  does  not 
represent  the  real  electrochemical  relationship,  as  it  is  described 
with  our  model,  leading  to  erroneous  results. 

4.2.  Saturation  and  hysteresis 

For  the  complete  model,  the  effect  of  the  saturation  concen¬ 
tration  X  and  the  hysteresis  remains  to  be  included.  To  this  end,  we 
consider  that  the  concentration  saturation  at  the  interface  affects 
the  charge  current  in  the  sense  that  when  X  is  equal  to  one/zero, 
only  a  fraction  of  the  charge  current  can  be  accumulated/delivered. 
In  such  circumstances  the  assumption  that  the  charge  current  is 
approx  equal  to  the  battery  terminals  is  not  fulfilled.  In  order  to 
obtain  such  current,  the  following  algorithm  is  used: 

if  X  >  1  and  fbat  >  0,  then  I  =  max(/bat  -  A(X  -  1),  0); 
else,  if  X  <  0  and  /bat  <  0,  then  I  =  min(/bat  -  AX,  0) 
else,  I  =  fbat;  end;  end; 

(45) 

where  /bat  is  the  battery  current  and  A  is  a  positive  constant  high 
value.  This  algorithm  takes  a  proportion  of  the  load  current,  so  that 
the  concentration  X  does  not  deviate  from  values  close  to  satura¬ 
tion.  Common  values  of  A  range  from  20  to  40. 

The  hysteresis  problem  arises  from  the  fact  that  the  EMF  in  the 
charge  is  different  from  that  obtained  in  the  discharge.  In  certain 
batteries  this  phenomenon  is  significant  and  must  be  taken  into 
consideration.  Assuming  f[j(0,X)  and  f^(O.  X)  are  the  charge  and 
discharge  of  EMF,  the  following  algorithm  was  used  in  our  model: 


Fig.  7.  Experiment  1.  Lithium  ion,  10  Ah. 
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fE(0,X(k))  =  p(k)fEc(0,X(k))  +  Cl-p(k))fEd(0,X(k)), 
p(k)  =  Ap(k-  1)  +  (1  -  A)sg(/bat(k)), 


(46) 

(47) 


r  _  1  ^bat(^)  ^bat(^  1) 


for|/bat(fc)  -  /bat(k  -  1)1  >  » 


(50) 


where  sg(Ibat(k))  is  one  if  the  sign  of  Ib3t(k)  is  positive,  and  zero  if  it  is 
negative.  The  forgetting  factor  A  is  within  the  interval  [0,1]  and 
weighs  the  velocity  of  change  between  both  EMF  functions. 

5.  Results 


where  N  is  the  number  of  current  steps  and  e  is  a  threshold.  Using 
this  equation,  the  calculated  internal  resistance  is  an  averaged 
value.  In  order  to  identify  the  parameters  do  and  m,  we  first  obtain  X 
from  the  voltage  error  and  the  inverse  of  the  EMF,  as  follows: 


In  this  section  we  present  the  identification  of  the  model  pa¬ 
rameters,  as  well  as  the  model  performance.  We  consider  the  sec¬ 
ond  order  approximation  of  equation  model  (36)  with  the  voltage 
model  given  by  the  approximation  in  equation  (35).  The  saturation 
and  hysteresis  are  also  considered.  The  model  parameters  to  be 
identified  are  Qmax,  do,  m,  and  req.  To  this  end,  the  continuous-time 
model  (31 )  is  discretized  using  a  sampling  period  h  =  60seg  and 
approximating  the  derivative  by  the  backward  shift  approximation, 
dx(t)/dt  ~  (x(kh)  -  x(kh  -  h))/h,  giving  the  following  discrete-time 
model: 

X(k  + 1)  =  0dx(k)  +  rldsoc(k)  +  r2dib3t{k)  +  X(k),  (48) 

where  $d  =  (1  -  hdo),  Tid  =  hd0,  and  r2d  =  /mj/CW.  The  dis¬ 
cretized  voltage  model  is  given  by 

EbatW  =  M0,X(k))  +  reqIbat(k).  (49) 

First,  the  value  of  req  is  obtained.  Taking  into  account  equation 
(48),  the  resistance  req  can  be  obtained  for  fast  variations  of  po¬ 
tential  due  to  increments  of  current,  as  follows: 


X(k)  =  fx( 0, Ebat(k)  -  reqIbat(k)) .  (51 ) 

In  the  case  where  the  EMF  has  hysteresis,  an  intermediate  value 
from  charge/discharge  values  was  considered  for  obtaining  X(k).  By 
assuming  that  the  value  of  Qmax  has  been  previously  identified  and 
that  the  disturbances  x(k)  are  small  with  respect  to  the  states,  using 
the  experimental  values  of  Jbat(k)  and  SO C(k)  together  with  the 
estimated  values  ofX(k),  the  coefficients  do  and  nj  can  be  identified 
by  using  equation  (48).  The  problem  leads  to  the  following  linear 
regression  model: 

X(k)-X(k-  1)  =  hd0(S(k-l)-X(k--l))+^lbat(k-l). 

(52) 

In  order  to  obtain  the  voltage/current  register  for  parameter 
identification,  we  propose  performing  two  tests.  On  one  side, 
starting  with  a  full-charge  battery,  a  train  of  regular  discharge 
pulses  of  current  is  used  to  determine  the  capacity  Qmax-  Also,  the 
SOC  is  obtained.  Assume,  as  well,  that  the  pulses  are  enough  spaced 
in  time  so  that  the  system  relaxes  between  pulses.  Taking  into 


Fig.  8.  Experiment  2.  Lithium  ion,  10  Ah. 


R.H.  Milocco  et  al.  /  Journal  of  Power  Sources  246  (2014)  609-620 


617 


Fig.  9.  Experiment  1.  NiMH,  Sanyo  2.7  Ah. 


Fig.  10.  Experiment  2.  NiMH,  Sanyo  2.7  Ah. 
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Fig.  11.  Experiment  1.  Lead  acid,  7.5  Ah. 


Fig.  12.  Experiment  2.  Lead  acid,  7.5  Ah. 
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Table  1 

Identified  parameters. 

Type  d0  [s_1]  re  [Ohms]  n,  [s_1]  Qmax  [Ah] 

Li  2.28  X  1(T3  1.13  x  1(T2  2.66  11 

NiMH  3.42  x  1CT3  2.05  x  10“'  5.67  2 

Pb  2.62  x  10“3  5.5  x  KT1  9.11  2.5 


account  that  after  the  system  relaxes,  the  value  of  X  is  equal  to  the 
SOC,  the  EMF  function  is  obtained  by  the  pair  (SOC,£)  measured  at 
the  end  of  the  relaxing  time  for  each  pulse.  In  order  to  identify  the 
hysteresis,  the  same  procedure  is  done  for  charge  pulses  to  obtain  a 
second  EMF. 

A  second  test  is  performed  by  using  a  random  sequence  of 
current  pulses  of  different  amplitudes  and  duration  in  both  charge 
and  discharge.  This  experiment  is  used  to  identify  the  constants  d0, 
and  ni  using  equation  (52). 

To  achieve  the  initial  state  of  charge,  SOC  =  0,  a  discharge  cur¬ 
rent  of  small  amplitude  is  applied  until  the  potential  decreases 
abruptly.  Since  the  current  is  small,  there  is  no  voltage  drop  across 
the  internal  resistance.  Consequently,  the  potential  variation  is 
solely  due  to  the  variation  of  surface  concentration  in  the  nonlinear 
region,  where  the  value  of  X  is  approximately  zero.  A  similar  pro¬ 
cedure  is  used  to  obtain  SOC  =  1. 

In  Figs.  7  and  8  both  experiments  are  shown  on  a  Lithium  ion, 
TM  SEIDEN,  Mode:  Prismatic  LiFeP04F  =  3.5  V,  nominal  capacity, 
Qmax  =  10  Ah  battery.  In  Figs.  9  and  10,  the  same  procedure  was 
used  for  a  Ni-MH,  TM  SANYO,  Model  HR-4UE  =  1.2  V  battery  of 
nominal  capacity  Q  =  2.7  Ah;  and  in  Figs.  11  and  12,  the  same 
procedure  was  applied  to  a  Lead  acid  battery,  TM  Press,  Model 
PR1240,  £  =  12  V  with  Q- 4  Ah  of  nominal  capacity.  The  parameter 
set  obtained  for  the  three  batteries  is  shown  in  Table  1.  Experiments 
have  been  performed  with  current  values  lower  than  1C.  Presum¬ 
ably,  the  resistance  values  do  not  vary  significantly  when  higher 
current  is  applied.  Should  this  not  occur,  a  new  parameter  identi¬ 
fication  at  high  currents  is  suggested. 

It  is  interesting  to  analyze  the  differences  of  the  RCE  among  the 
three  batteries  tested.  Solving  equation  (36)  for  a  constant 
discharge  current  entry,  I  —  -mQmax.  where  m  is  an  integer  in  units 
(1/time),  the  states  SOC(t)  andX(t)  can  be  obtained.  By  performing 
the  subtraction  of  the  two  states,  the  loss  of  capacity  due  to  the  RCE 
is  obtained  as 


Fig.  13.  Rate  capacity  for  the  three  types  of  batteries. 


SOC(t)  - X(t)  =  m"1^  (l  -  e~dot) .  (53) 

It  can  be  seen  from  the  equation  that  the  parameters  do  and  nj  of 
the  HPF  define  the  RCE.  In  Fig.  13  the  capacity  lost  due  to  RCE  of  all 
the  batteries  tested  is  depicted  for  constant  current  with  m=  1  h 
It  can  be  seen  that  the  Lithium  battery  has  less  RCE  than  Ni-MH,  and 
it  can  also  be  seen  that  the  lead  acid  battery  shows  a  great  RCE 
effect. 

6.  Conclusions 

A  simple  model  based  on  an  electrochemical  approach  to  be 
used  in  applications  BMS  has  been  presented.  The  key  feature  is 
that  it  can  be  used  for  any  type  of  battery  that  works  for  charge 
accumulation  and  transport  of  matter.  The  model  is  based  on  two 
fundamental  electrochemical  states  as  the  state  of  charge  and 
the  surface  concentration.  The  latter  is  responsible  for  the  RCE 
and  enables  a  clear  understanding  of  potential  variation  as  a 
function  of  the  surface  concentration  X.  Tests  show  that  a  low-order 
approximation  is  indeed  sufficient  to  obtain  a  good  model. 
Furthermore,  we  have  analyzed  the  similarities  and  differences 
with  respect  to  other  widely  used  simple  models,  demonstrating 
that  our  model  has  the  same  degree  of  simplicity  associated  with 
important  electrochemical  variables,  which  expands  its  range  of 
use. 
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Appendix  A 

The  inverse  function  £  — /e(/,X)  of  equation  (5)  exists.  In  order  to 
show  this,  let  us  state  that  a  continuous  function /is  one  to  one  if 
and  only  if  it  is  either  strictly  increasing  or  decreasing  with  no  local 
maxima  or  minima.  In  our  case,  the  partial  derivative  with  respect 
to  £  after  some  algebraical  manipulation  can  be  written  as 

0X  _  k°ba0ea°E(l  -X)  -  krarea'EX 

dl~  krea'E  +  k°ea°E  '  1  } 

Taking  into  account  that  X  e  (0,1 ),  and  ar  <  0,  a°  >  0,  the  de¬ 
rivative  above  is  always  positive.  Thus,  for  each  value  of  /,  the  in¬ 
verse  function  £  =  /e(X,7)  exists.  Hence,  there  is  only  one  value  of  £ 
for  each  couple  (X,/). 

Appendix  B 

By  taking  into  account  that  S  =  a  +jto  the  infinite  dimensional 
rational  function  H(S)  can  be  written  as 

H(c r  +jv)  =  BJ((oJ  -  A)  +j(iil)  1  yBF.  (55) 

Thus,  when  «  =  0,  H(<r)  =  BT  (ol  -  Aj^'yBF,  which  is  real  for  all 
(j,  the  first  condition  is  fulfilled.  To  prove  the  second  condition, 
taking  into  account  that  f  and  y  are  positive  real  constants,  we  find 
that  H(o  >  0)  >  0  iff  -A  >  0.  Thus,  we  need  to  show  that  -A  is 
positive  semidefinite,  which  is  the  same  as  proving  that  given  an 
arbitrary  raw  vector  Z  always  fulfills  that  —ZAZT  >  0.  Note  that  each 
element  ky  appears  in  matrix  A  at  positions  ( ij )  and  Qi)  with  a 
positive  sign,  and  at  (jj)  and  (ii)  with  a  negative  sign.  Then,  each 
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element  ky  appears  multiplied  by  — (Zj  -  Zj)2,  which  is  always 
negative  semidefinite  for  arbitrary  values  of  Zj  and  Zj.  Since  the 
elements  ky  are  positive,  we  conclude  that  the  matrix  -A  is  positive 
semidefinite  and  the  second  condition  for  PRF  of  H(S)  fulfills. 

Appendix  C 

Assume  the  potential  E  of  the  Butler-Volmer  equation  in 
equation  (4)  is  a  sum  of  two  components.  One  due  to  the  EMF  given 
by  equation  (6),  which  does  not  depend  on  the  current,  denoted  as 
Ei  — /e(0,X),  and  the  other  which  depends  on  the  current  I,  called  Ej. 
The  functional  relationship  between  E2  and  the  current  can  be 
obtained  by  replacing  E  in  equation  (4)  which  gives 

I  =  Y(e“°£2  -e“r£2),  (56) 

where  y  =  k°Xea"E'  =  kr(l  -X)e“r£>  come  from  the  fact  that  E,  is 
defined  for  /  =  0  and  £2  =  0.  Then,  the  voltage  equation  (34)  can 
now  be  written  as  follows: 

4at  =/E(0,X)  +  (£2  +  re/bat),  (57) 

The  second  term  on  the  right  can  be  interpreted  as  the  ohmic 
drop  produced  by  an  equivalent  resistance  req(y,E2)  given  by 

req(y,E2)  =  ^-+re  (58) 

Now,  we  want  to  know  the  dependency  of  re q(.y,E2)  with  both 
current  and  concentration.  To  this  end,  note  that  the  symmetry 
factor  of  electrodes  in  a  rechargeable  battery  has  a  value  approxi¬ 
mately  equal  to  a0  =  0.5.  This  means  that  the  charge  current  cor¬ 
responding  to  positive  overpotentials  are  approximately  the  same 
than  for  discharge  current  corresponding  to  negative  over¬ 
potentials.  Then,  the  equality  a0  =  —ar  holds  in  equation  (56)  and 
the  current  can  be  written  as 

l  =  2y  sinh(a°E2),  (59) 

Thus,  by  assuming  /bat  =  /,  the  nonlinear  resistance  is  given  by 
E2  1  sinh-1  (J/2y) 

I  2ya°  I/2y  '  {  * 

It  must  first  be  noted  that  the  element  sinh-1  (x)/x  is  a  sym¬ 
metrical  decreasing  function  with  a  maximum  value  equal  to  one, 
atx  =  0.  Thus,  the  nonlinear  resistance  has  it  maximum  value  equal 


to  1 12ya°  at  /  =  0.  As  the  current  increases  the  resistance  decreases. 
It  is  also  worth  mentioning  that  the  constant  y  increases  with 
battery  capacity  as  may  be  inferred  from  equation  (56).  Then,  it  is 
expected  that  for  large  capacity  batteries  the  non-linear  distortion 
is  less  significant.  Finally,  it  is  easy  to  see  from  equation  (60)  that  for 
any  given  value  of  current  amplitude,  the  nonlinear  resistance  de¬ 
creases  as  y-  or  equivalently  X-increases. 
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